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O": ABSTRACT 

o . 
o ■ 

We study the stochastic coagulation equation using simplified models and 
Q_i. efficient Monte Carlo simulations. It is known that (i) runaway growth occurs 

if the two-body coalescence kernel rises faster than linearly in the mass of the 
heavier particle; and (ii) for such kernels, runaway is instantaneous in the limit 
! that the number of particles tends to infinity at fixed collision time per particle. 

Superlinear kernels arise in astrophysical systems where gravitational focusing is 
! important, such as the coalescence of planetesimals to form planets or of stars 

to form supermassive black holes. We find that the time required for runaway 
decreases as a power of the logarithm of the the initial number of particles. 
^ Astrophysical implications are briefly discussed. 

! Key Words: planetesimals, planetary formation; collisional physics; methods, 

numerical. 
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1. Introduction 

> ■ 
• i— i . 

A frequently-encountered process in many fields of science is the random coalescence 
of small bodies into larger ones, conserving total mass. Astrophysical examples include the 
coalescence of planetesimals into planets (Safronov 1969) and of stars into black holes (Lee 
1987; Quinlan and Shapiro 1990). When the number of bodies is large, coalescence is often 
modeled by Smoluchowski's equation (Smoluchowski 1916), also known as the statistical 
equation: 

dc 1 °° 

-7T = 2 K ( l ~ ■'■■' ]r ' ■>''■> ~~ C *J2 K ( i ^)c j , (1) 

T 3=1 3=1 

Here the concentration Cj(r) is the number of bodies per unit volume of mass m,; oc i at 
time r, and the functional form of the coagulation kernel K(i,j) = K(j,i) is chosen to 
approximate the mass dependence of the two-body collision rate. 



The statistical equation is intended to model coalescence on average, smoothing over 
fluctuations. It is not expected to be accurate when the initial number of particles, N, is 
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small. More interestingly, eq. (1) can fail even for iV 3> 1 if K(i,j) increases sufficiently 
rapidly with % and j. Coalescence should conserve mass, so that the mass density 

oo 

p(j) = J2kc k (r) 
k=l 

is constant. It is easy to show from eq. (1) that p = provided that the relevant summations 
converge and can be interchanged. But in the analytic solution of the case K(i,j) = ij 
starting from the monodisperse initial conditions 

N 

c i (0)=6 il , p(0) = - = l, (2) 

p(r) begins to decrease after r = 1 (Trubnikov 1971). This is usually interpreted to mean that 
a macroscopic runaway particle has formed — also known as a gel, because of applications 
in physical chemistry. It is believed that for K(i,j) ~ (ij) v with v > 1, gelation begins 
immediately in the statistical equation (Jeon 1999; Lee 2000, and references therein). The 
model (1) can be extended to include the gel/runaway particle explicitly (Flory 1953). Monte 
Carlo simulations for large but finite N show, however, that such models do not accurately 
predict the time dependence of the gel mass even on average if v > 1 (Spouge 1985). 

Coagulation is more correctly described by the stochastic equation for the joint proba- 
bility /(rii, n 2 , ■ ■ ■ ; t) for the occupation numbers {rii, n 2 , . . .} of the mass bins: 



d N N 

f( n ',t) = Y Y K(iJ)(n i + l)(n j + l)f(...,n i + l,...,n j + l,...,n i+j -l,...]t) 

i=l j=i+l 



dt 



1=1 

N N N ^ 

i=l j=i+l i=l 

The statistical equation results from taking first moments, 

Ci(t) = V~ l ni(t) = V' 1 J riif(ni, n^i, rii, n i+1 , . . . ; t) dnx . . . dn i - 1 dn i+1 (4) 

in the limit that N, V — > oo at fixed p(0). In order to get a closed set of equations, one 
assumes that rQfTj = UiX nj. This is justified if the occupation numbers are approximately 
uncorrelated. But if runaway should occur, then the occupation numbers of all the low-mass 
bins are correlated with bins i ~ N traversed by the runaway particle. 
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Like the statistical equation, the stochastic equation has analytic solutions for the kernels 
K(i,j) oc constant, % + j, and i x j (Lushnikov 1977; Tanaka and Kiyoshi 1993). For the 
first two kernels, the predictions of the two equations are compatible in the following sense: 
after monodisperse initial conditions, q(t) computed from eq. (1) agrees with computed 
from eqs. (3)- (4) in the limit iV — > oo if t and r are related by the initial collision time per 
particle, 

toon = [(N - 1) K(l, I)]' 1 , r = t/t coll . (5) 

For these kernels, there is no runaway, and the statistical equation conserves mass. Even for 
K(i,j) = ij, a similar agreement is found between the statistical and stochastic equations at 
t < 1 and mass bins i <C y/N (Tanaka and Kiyoshi 1993). The stochastic equation confirms 
that runaway begins at r > 1; more precisely, Lushnikov (1977) shows that the quantity 

N 

»=i 

which can be interpreted as the runaway mass fraction on average, becomes nonzero at r = 1 
and pal — 2e~ T at r >■ 1. 

Kernels such that K(i,j) oc % v when % ^> j are sometimes called "unphysical" if v > 1 
because if collision rate were proportional to surface area, then surface area would have to 
increase more rapidly than mass. But in astrophysics (and elsewhere), gravitational focusing 
and mass stratification may conspire to produce v > 1 (§5). Thus a better term for such 
kernels would be superlinear. The occurrence and timescale of runaway are of great interest 
since the runaway particle may represent a planet or supermassive black hole, for example. 
Unfortunately, few general results are known for the superlinear regime. Spouge (1985) 
conjectured on the basis of Monte Carlo simulations of stochastic coagulation that runaway 
is instantaneous for v > 1 in the limit N — > oo: in other words, the entire mass is consumed 
by a single particle after an infinitesimal multiple of the single-particle collision time t co ii. A 
similar conjecture had previously been made by Domilovskii, Lushikov, and Piskunov (1978) 
for the special case v — 3, also on the basis of Monte-Carlo simulations. Recently, Jeon 
(1999) has supplied a proof. 

Clearly, runaway/gelation is not instantaneous for finite N. To the best of our knowl- 
edge, there have been no general and quantitative statements about the scaling of the run- 
away time with N, although the question is an important one since N is never truly infinite 
in practical applications. The main result of the present paper will be the conjecture, sup- 
ported by Monte Carlo simulations, that r runaway varies as a negative power of logiV for 
v > 1. Even before performing our simulations, we were lead to this conjecture on the basis 
of the highly simplified model presented in §2. Our Monte Carlo algorithm for stochastic 



k(t) == lim 
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coagulation is described in §3 and tested against the analytic solutions of the stochastic 
equation cited above. Simulations for v > 1 are reported in §4. Finally, in §5, we discuss 
applications of our conjecture to cases of astrophysical importance. 



2. MONOTROPHIC MODEL 

The full stochastic equation seems to be analytically intractable for general values of 
the merging exponent v. In this section we solve a simplified problem in which a single 
"predator" feeds upon an unevolving population of unit-mass "prey." The prey do not 
merge with one another but only with the predator, whose initial mass is equal to that of 
the prey. A slight additional simplification results from assuming that the prey population 
is infinite. We call this model "monotrophic" after the Greek mono (one) + trophein (to 
nourish). Our Monte-Carlo simulations of the full stochastic equation are well described by 
the monotrophic model in their later stages when (or if) the most massive particle exhibits 
runaway growth. 

We work in time units r such that the feeding rate of a unit-mass predator is unity, and 
we assume that the feeding rate increases as the z/ th power of its mass. Let the predator start 
with mass equal one at time r = 0, and let p(k, r) be the probability that it is in the k th mass 
bin and correspondingly has mass = k after time r. The probability that the predator will 
enter mass bin k during the brief time interval (r, r + dr) is clearly p{k — 1, t) x (A; — 1)" x dr, 
and the probability that it will vacate bin k during the same interval is p{k, r) x k v x dr. 
Hence the evolutionary equation and initial conditions for p are 

^-(k,T) = (k-iy P (k-l,T)-k v p(k,r), and p(k,0) = 6 ktl . (6) 

d °° 

One might suppose that — > p(k, r) = 0, since the sum of the righthand sides (6) would 

(IT < J 



k=i 



appear to cancel, but after solving eqs. (6), one finds that 

oo 

*(r)^l t) (7) 
k=i 



can be nonzero for r > 0. Clearly 0(r) is the probability that the predator is not to be 
found in any finite mass bin. So we interpret <f>(r) as the probability that runaway growth 
has occurred. From eq. (6), 0(0) = and 

^(r) = lim k v p(k,r). (8) 

aT fc^oo 



- 5 - 



Hence, the runaway cannot occur until the mass probability function develops a power-law 
tail p(k, r) oc k~ v . But Y^k=iP(k, t) must be finite (in fact < 1), and the power-law tail has 
a convergent sum only if v > 1. So we conclude that runaway occurs in the monotrophic 
model only if v > 1. 

To solve for 0(t), we take the Laplace transform r — > z of eq. (6), 

zp(/c, - 4,i = (k - l) u p(k - 1, z) - k u p(k, z). 

This is easily solved: 

k _ 1 

fc"^^) = n (i + (*>i). (9) 

In view of eq. (8), the Laplace transform of (f) is 

^) = ;£5» n( 1 + F)~'- (10) 

A subscript has been placed on <j> v (z) to acknowledge its dependence on v. Standard con- 
vergence tests show that the limit of the product vanishes if v < 1 except at the poles 
z = —l u , —2 U , —3^, . . ., so that 4> u <i(t) = also vanishes. Once again, therefore, we see that 
runaway does not occur for v < 1. On the other hand, if v > 1 then the product has a finite 
nonzero limit. 

Exact results are possible for v = 2. The product (10) has the closed form 

Mz) = . I (n) 

y/z smh(ir^/z) 

The only singularities of this function at finite z are simple poles on the negative real axis 
where \f—z is a positive integer. Evaluating the inverse Laplace transform by residues yields 
an infinite series that converges rapidly at large r: 1 

oo 

2 (r) = l + 2^(-l)V^ (12) 
k=i 

Writing 

oo 

2 (r) = J A{x)e~ Tx2 dx, 

— OO 

OO oo 

where A(x) = ^ (-l) k 5{x - k) = 2 cos [(2n + l)irx] , 

fc=— oo n=0 



1 4>2{t) can also be expressed in a parametric closed form involving elliptic integrals. 
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and integrating the Fourier series term by term yields a series that converges rapidly at small 
times: 

7r 2 (2n + l) 2 " 



/ — oo 

<h( T ) = 2 J~ $^ exp 

' T n=0 



At 



(13) 



Eq. (13) shows that runaway is exponentially unlikely for r 1, while eq. (12) shows that 
it is virtually certain for r 3> 1. The probability reaches 50% at r « 1.37. 

Similar statements hold for general values of v > 1 (cf. the Appendix): 

(27r)^- 1 )/ 2 v /^7F = !) r ex P {~iy ~ 1) [^sin(7r/z/)/7r]^ /( ^ 1) r" 1 /^- 1 )} if r < 1, 
&(r)«< ' (14) 



1 - 



ffi°= 2 (l-^) exp(-r) ifr»l. 



These expressions agree with the leading terms of the series (12) and (13) when v — 2. 
Scrutiny of the first exponential in eq. (14) reveals that as v approaches unity from above, 
the time at which the runaway probability first becomes appreciably different from zero is 
~ (u- 

In the original model of §1, all particles can be predator or prey. For r « 1, the 
results of this section suggest that the probability that any given particle has achieved large 
mass is ~ 4 > { T )- We presume these probabilities to be approximately independent if N is 
sufficiently large, because two randomly selected particles could both grow to a mass ^> 1 
before competing for the same prey. Thus the total probability of runaway at a time r << 1 
is pa N(j){r). Given the estimate (14) for 0(r), the total runaway probability approaches 
unity at a time scaling as 

T„a y oc(logiV) 7M (15) 

with 7(1/) = 1 — v < 0, but only if v > 1. To test this conjecture, we undertook the Monte 
Carlo simulations described below. 



3. Numerical Technique and Tests 

Direct numerical solutions of the stochastic coagulation equation (3) are difficult to 
obtain for large N, so we have used Monte Carlo simulations. Great savings in memory and 
processor time can be achieved by keeping track of the occupied bins only. Therefore, we 
adopt an indexing scheme different from that of §1: at each time step, rii 7^ will represent 
the occupation number of the i th nonempty bin, with mass per particle and 1 < i < I. 
For expository convenience, the bins are sorted by mass (rrii < rrij if i < j), although this is 
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not essential to the algorithm. The total mass is 

^n imi = N. (16) 



i=i 



The storage required by our scheme is 0(1). Because mj > i and rij > 1, the summation 

1 /2 

(16) is > /(/ + l)/2. Hence the number of occupied bins / < (27V) 1 . In our simulations, 
the maximum / encountered is almost always much smaller than this upper bound. 

The coagulation kernel is now written as K(mi,mj) rather than K(i,j) to emphasize 
that the coalescence rate depends upon the particle masses rather than the arbitrary bin 
indices. The total coalescence rate of bins i, j > i is 

def ,., _ N f I'," j i < j 



RiJ = K(m i ,m j )x^ ni{ni _ 1)/2 . = j . (17) 

It is useful to imagine these rates arranged as an / x / upper-triangular matrix, and to 
introduce the partial sums 

i-l / j 

= ^2^2 R k,e + ^2 Ri,e, (18) 
fc=i i=k e=i 

i.e. Sij is the sum of the first i — 1 rows plus the first j columns of row i, and the total 

def 

coalescence rate S = Sij. We store in the computer memory the rates Sij, % = 1,2,. ..I. 
Each merging alters the occupation numbers of at most three bins, so these rates can be 
updated in 0(1) operations. Indirect indexing minimizes the work of inserting or deleting 
rows and columns as the list of occupied bins changes. 

Every simulation begins with monodisperse initial conditions, 

m = N, mi = 1, 1 = 1, (19) 

and ends after N — 1 mergings with 

ni = 1, m 1 = N, I = 1. (20) 

Merging occurs at random but increasing times ti < ti . . . < tN-i chosen as follows: If S is 
the total coalescence rate computed after the s th merging, the probability that no further 
merging occurs before t > t s is 

V = exp [-S ■ (t - t s )} . 

Therefore, we chose a random number X ~ U(0, 1] (i.e., X is uniformly distributed between 
and 1) and take 

*,+i = * a -S'HnX. 
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Thus merging occurs at t 8+ i, and the next task is to decide which bins are involved. Choosing 
a second random number Y ~ [7(0, 1], we find i, j (j > i) such that 

Sij-! < Y-S < Sij. (21) 

Since S it j — Sij-i = R i: j [eq. (18)], bins i,j are selected with the correct probability Ri^/S, 
all rates being evaluated just before the s + 1 st merging. 

With this scheme, there are N — 1 steps per simulation and 0(1) operations per step, 
so the number of operations per simulation is at most 0(N 3 ^ 2 ). In practice, / <C N 1 ^ 2 so 
that the computer time is almost linearly proportional to N. We average many simulations 
to obtain the statistics of quantities of interest. 

We have tested our code against the exact solutions of Tanaka and Kiyoshi (1993) for 
the cases K(m^m^) = 1 and K(m^m^) = rrii + m,j. We compare at N — 10 and find 
excellent agreement (Fig. 1). Not surprisingly, even though runaway does not occur, the 
statistical model is a poor approximation at such small N. 



4. SIMULATION RESULTS 

We discuss simulations for two classes of coalescence rates, K(rrii,mj) = (rrii x rrij) v 
(multiplicative kernel) and K(mi,rrij) = (rrii + mj) u , v > 1 (additive kernel). 

Figure 2 shows the evolution of various averaged diagnostics for the representative case 
v — 2, N = 10 7 and for both additive and multiplicative kernels. Instead of elapsed time, we 
use the logarithm of the number of mergings, A^ merg , as the independent variable. Fig. 2(a) 
shows the (average) relationship between these coordinates. Note that we have reintroduced 
t [cf. eq. (5)], which is time normalized to the collision time per particle in the initial state 
(19). We see that there are two evolutionary phases. In the first phase, the total coalescence 
rate is dominated by mergings among unit-mass particles [Fig. 2(d)]. The mass spectrum, 
i.e. the distribution of average occupation number n with particle mass m, is continuous, as 
shown by the near-equality of the largest and second-largest mass [Fig. 2(b)], and extends 
to steadily higher m [Figs. 2(e,c)]. The total coalescence rate is approximately equal to its 
initial value N/2t co \\. This is the statistical growth phase. 

The runaway growth phase could be considered to begin when there is a significant 
mass gap between the highest and second- highest occupied bin [Fig. 2(b)]. We find it more 
convenient to declare runaway when the rate of mergings involving the highest occupied bin 
is 50% of the total coalescence rate [Fig. 2(d)]. At this point mj, the mass per particle of 
the highest-mass bin, is as N l / U . The runaway phase accounts for a negligible fraction of the 
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total time (t or r) but a majority of the total mergings [Fig. 2(a)]. 

For other coalescence exponents v > 1, the qualitative features of the evolution are 
similar to those displayed in Fig. 2 at sufficiently large N. As v approaches unity, larger N 
is required to obtain similar behavior, and runaway begins later — i.e., at larger N meige /N. 

At a given N and u, runaway occurs later for additive than for multiplicative kernels and 
involves more particles of intermediate mass. The latter is explained by the fact that the col- 
lision rate between particles of very different mass is almost independent of the lighter mass 
with the additive kernel. Once the runaway phase of the multiplicative case is well estab- 
lished, most of the particles are of unit mass, and the evolution is similar to the monotrophic 
model [Fig. 2(e,f)]. 

Figure (3) shows that the runaway time slowly decreases with increasing N, even after 
scaling by the single-particle collision time (5). At sufficiently large N, it appears that 
^runaway decreases as a power of log N and that the slope is steeper for larger u, in qualitative 
agreement with the conjecture (15). The first panel of Fig. (4) shows rough quantitative 
agreement with the conjecture in the slope of "f(v) against v for multiplicative kernels. 
There appears to be an offset in the vertical intercept, but one knows from the exact results 
of Lushnikov (1977) that r runaway is independent of N as N — > 00 at v — 1, so the solid and 
dashed lines should intersect at 7(1) = 0. Unfortunately, it is extremely difficult to measure 
7 near v — 1 from our simulations. The second panel of Fig. (4) shows that for additive 
kernels, the empirical slope of 7(1/) is close to one half of the prediction (15). 

5. DISCUSSION 

Our Monte Carlo simulations indicate that runaway is instantaneous in the limit N — > 00 
if the coalescence kernel is superlinear, in agreement with previous authors (Jeon 1999; 
Spouge 1985). But the simulations, as well as the heuristic model of §2, suggest a more 
quantitative new result: namely, that the limit r mnaway — > is approached only logarithmi- 
cally in N for power-law kernels. 

We consider briefly two astrophysical applications: coalescence of planetesimals into 
planets in a protostellar disk, and black hole formation through stellar mergers in galactic 
nuclei. In both cases, the two-body collision rate per unit volume is ^^(av)^, where z/j is 
the number of particles of mass m ; per unit volume, and the rate coefficient is approximately 
(Lee 2000) 

(crv)ij = Tx(Ri + Rj) 2 



Vij + 



2G{mi 



rrij 



(Ri + Rj) 



■'i.i 



(22) 
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Here Ri is the radius corresponding to mass m^, and % is the root-mean-square relative 
velocity of the particles. The prospects for runaway are usually more favorable when the 
second (gravitational focusing) term in the square brackets above is dominant. Henceforth, 
we neglect the first term. A power-law mass-radius relation 

R oc m a (23) 

describes planetesimals if a ~ 1/3, and stars if a ~ 1. For a maxwellian velocity distribution 
with equipartition of kinetic energies among mass groups, 

I + J X 



Vij = Vi 

\ ij 

where v i is the velocity dispersion of particles of mass m x and mj = i x mi. Then 

H« - (*" + r) (< + j) 1/2 (^) 1/2 (24) 

Since this is proportional to j a+1 for j ^> i, runaway is possible if a > 0. 

Mass stratification somewhat intensifies the tendency towards runaway. A large-scale 
gravitational potential will concentrate the heavier particles into a smaller volume than the 
lighter ones. Thus in a disk of planetesimals with orbital angular velocity Q, the distribution 
of mass group i with height z above the midplane is 

t>i(z) w i/i(0) exp (-z 2 /2H?) , ^ = ^/fi oc r 1/2 . 

Hence upon integrating with respect to z, 

1/2 



/ 



Here hi is the number of particles of mass m, per unit area. The mass-dependent factor 
in this expression multiplies the rate coefficient (24) to yield an effective collision kernel 
(a w 1/3) 

K(i,j)=K(l,l) x ijpfi+j 1 ' 3 )^, (26) 

which has the same j 4//3 scaling for j 3> % as the rate coefficient (24) but a stronger depen- 
dence on the mass of the lighter particle: % versus i 1 ^ 2 . This is likely to cause more rapid 
consumption of intermediate-mass objects, so we expect evolution as if for a multiplicative 
kernel [K(i,j) oc (ij) 4 ^ 3 } rather than an additive one [K(i,j) oc (i + j)^ 3 ]. 
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The situation in galactic nuclei is more complicated if the heavier stars concentrate into 
a selfgravitating population. Setting aside this possibility (which would make runaway still 
more rapid) and assuming a harmonic mean potential, we get in place of eq. (25) 

o 

since the geometry is spherical rather than planar. Here is the total number of particles 
of type i. So in this case, with a — 1, 

K(i,j) = K(l,l) x (ij) 2 . (27) 

This is precisely our multiplicative v — 2 kernel. 

As many as ~ 10 11 asteroid-sized planetesimals are required to make up the mass of a 
terrestrial planet. The smallest coalescence exponent for which our simulations give reliable 
averages for r runaway is v — 3/2, and the largest N we have used is 10 9 . Nevertheless, by 
extrapolating the curves in Fig. (3), we estimate t runaway ~ 0.4t coll at N = 10 11 , versus 
Runaway ~ 2 i co ii at TV = 10 2 . The dependence on logiV is somewhat stronger than one would 
expect from the asymptotic scaling (15) because of the curvature in Fig. (3). 

Black holes associated with quasars in the nuclei of galaxies are thought to have masses 
of order 10 8_9 M Q . The prevailing wisdom is that the black hole forms by accretion from 
a gas disk, but formation by stellar coalescence has long been a popular alternative model 
(Rees 1984, and references therein). Referring again to Fig. (3), but now for v = 2, we find 
Runaway ~ .09 t co n at N = 10 9 versus trunaway ~ .35t co ii at iV = 10 2 , roughly in agreement 
with the predicted asymptotic scaling (logiV) 1 ^. 
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A. Asymptotics of the monotrophic runaway probability 

To derive eqs. (14), one must approximate the inverse Laplace transform 

C+ioo 



W = 2^ / Mz)e ZT dT, (Al) 



C—ioo 



with (f> v (z) represented by the infinite product (10). The real constant C must be positive but 
is otherwise arbitrary, since the integrand is analytic for Re(z) > 0. Clearly, (j> v (z) decreases 
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rapidly as Re(z) > increases, while e 2T has the opposite trend. A good choice for C is the 
point C along the real axis where the product of these two factors in the integrand (Al) 
is smallest. By analyticity, this is actually a saddle point with respect to complex z, so the 
integrand will further decrease along the integration contour with increasing |Im(z)|. For 
T <C 1, Co ^> 1, because the exponential varies only slowly; whereas for r ^> 1, Co lies close 
to the pole of (f> v {z) at z — 0. 



Thus for small r, we are concerned with Re(z) ^> 1. We rewrite eq. (10) as 

k 



In 



zcj)(z) 



-5> 

fc=i 



lim 

n^oo 



k=l 



EM1 + T + ln 



{n\y 



(A2) 



The latter sum is dominated by its upper limit and can be approximated by the leading 
terms of the Euler-Maclaurin formula (Bender and Orszag 1978): 



5>(i + 



k=l 



In 1 



ar 
z 



dx + -In 1 + 



= in + - In 1 + — 



fi- 



nis + 1/2; 



dx 



z + x h 



where integration by parts has been used. As n — > oo, the latter integral becomes 



vz 



dx 



z + x u 



vz 



1/1/ p 



r l- 



7T 



sin(7r/z/) 



After this is substituted into eq. (A2) and Stirling's approximation is used for n\, the terms 
growing with n cancel, and upon taking n — > oo one has 

ln0(2) w - , — - -lnz + -ln(27r) +0(^ 1 ). (A3) 



Using this, we estimate 



Co 



"/("-!) 



rz/ sin(7r/f ) 

and then perform a straightforward steepest-descent contour integration, i.e. we expand 
In 0(2) + zt as a Taylor series to second order in z about z = C and use the resulting 
Gaussian to approximate the integrand (Al). This gives the upper line in eq. (14). 

For r > 1, we evaluate (Al) from the residues of the poles at z = and z = —1; the 
remaining residues are strongly suppressed by the exponential factor. This yields the lower 
line in eq. (14). 
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Fig. 1. — The number of particles, averaged over 1000 simulation runs and normalized to 
the initial number of particles (N = 10), versus our dimensional time is shown by the two 
solid lines. The analytical results for the stochastic and statistical coagulation equations are 
shown by the two dotted lines (which almost coincide with the corresponding solid lines) 
and two dashed lines respectively. 
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Average evolution for 300 simulations with N = 10 7 . Solid lines: K(m i: rrij) = 
(rrii x rrij) 2 ; dotted lines: K{m^m^) = (rrii + rrij) 2 . (a): time versus number of mergings, 
iV merg ; (b): the mass per particle of the highest (monotonic curves) and the second highest 
bins, i.e. rrij and mj_i; (c): the number of particles in the highest bin; (d): fraction of the 
total coalescence rate that involves one particle taken from the highest bin (monotonically 
increasing), and by mergings between unit-mass particles (decreasing); (e): number of occu- 
pied mass bins, J; (f): number of unit-mass particles as a fraction of the total remaining, i.e. 

ni/(N - Emerge)- 
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K(m i ,m j ) = (rri; x m J ) l/ K(m i ,m J ) = (rrij + m } ) 
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Fig. 3. — Average normalized time at which runaway begins versus initial particle number 
iV for different coalescence exponents v. Error bars enclose 50% of the simulations [30 — 10 5 
simulations at each (TV, v)\. 
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Fig. 4. — Exponent 7 for the relation r runaway oc (logiV) 7 between runaway time and particle 
number, versus coalescence exponent v. Solid lines: best fit based on asymptotic slope 
seen in Fig. 3 at large N (the slope and error bars are obtained by the least-squares linear 
regression). Dashed lines: prediction 7 = 1 — v of monotrophic model [eq. (15)]. 



